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1. Introduction 

The fractional calculus has allowed the operations of integration and differentiation to any 
fractional order. The order may take on any real or imaginary value. This fact enables us to 
consider the order of the fractional integrals and derivatives to be a function of time or some other 
variable. Lorenzo and Hartley [ 1 998 J first suggested the concept of variable order (or variable 
structure) operators. In this concept the order of the operator is allowed to vary either as a 
function of the independent variable of integration or differentiation (?) or as a function of some 
other (perhaps spatial) variable (v). In that paper a preliminary study was done on several 
potential variable order definitions where the order was a function of the independent variable of 
integration or differentiation (/), and initial properties were forwarded. 

This paper explores more deeply the concept of variable order integration and differentiation 
and seeks to create meaningful definitions for variable order integration and differentiation and to 
relate the mathematical concepts to physical processes. Several candidate definitions are 
developed and characterized. For these definitions operational methods are provided, and 
example applications are shown. One of the important characteristics of these operators is the 
manner in which the operator accounts for the order history (or order memory). Two types of 
memory are attributed to the operators, first w, , is the fading memory of the overall operator, and 
second is the memory related to the history of the order of the operator. This second type of 
memory, called nu , is found to be an important attribute of the operator. Measures of these types 
of operator memory are presented. 

From these candidate characterizations the concept of “tailored” variable order integration and 
differentiation is evolved. By this we mean the ability to choose the form of the operator 
definition to control how the operator deals with its own order history. 

The concept of the (continuous) order distribution, where the order is distributed over 
variables other than the independent variable was suggested by Hartley and Lorenzo [2000] and 
an initial and important application of order distributions was made to the identification of 
fractional order systems. In a related work Bagley and Torvik [2000] develop a special transform 
for the solution of distributed order differential equations. The present paper looks at two forms 
of order distributions and their applications to dynamic processes. Also, a variety of physical 
problems that motivate interest in these areas are briefly presented. Finally, the paper defines 
operators of multi-variable order and operators that combine the variable order and the order 
distribution concepts. 


*Distinguished Research Associate 
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2. Motivation for Variable Fractional Order Operators 

There is considerable potential physical motivation toward the creation and implementation of 
the concept of variable order operators. A few possibilities will be mentioned here. From the field 
of viscoelasticity (Bland 1960), the effect of temperature on the small amplitude creep behavior 
(force/extension) of certain materials is to change the characteristics from elastic (spring-like, 
q = 0 ) to viscoelastic or viscous (damper-like, q = -1 ). This relates to the expression 

C D;« F(r)=k(x a (t)-.x b ( r )\ (2.1) 

where F is force and x is displacement. Experience (and experiments), typically, is based on 
fixed temperature, but real applications may well require a time varying temperature to be 
analyzed. Polymer linear viscoelastic stress relaxation was studied by Bagley (1991). This 
process is described by fractional differential equations of order p for a given fixed temperature. 
The paper shows a clear dependence of p on temperature for polyisobutlene and correlates the 
fractional model and experiment. Further, it is indicated that p order fractal time processes lead 
to p order fractional derivative constitutive laws. Smit, W. and deVries (1970) studied the stress- 
strain behavior of viscoelastic materials (textile fibers) with fractional order differential equations 
of order a , with 1 > a > 0. They show, based on related experiments, a to be dependent on strain 
level. 

Glockle and Nonnenmacher (1995) studied the relaxation processes and reaction kinetics of 
proteins that are described by fractional differentia] equations of order p . The order was found to 
have a temperature dependence. 

Electroviscous or electrorheological fluids (Klass and Martinek) and polymer gels (Shiga) are 
known to change their properties in response to changes in imposed electric field strength. The 
properties of magnetorheological elastomers respond to magnetic field strength (Davis). 

From the field of damage modeling, it is noted that as the damage accumulates (with time) in a 
structure the nonlinear stress / strain behavior changes. It may be that this is better described with 
variable order calculus. Finally, the behavior of some diffusion processes in response to 
temperature changes may be better described using variable order elements rather than time 
varying coefficients. 

3. General Preliminaries for Variable Order Operators 

3.1 Issues 

This section presents several issues and characteristics that generally apply to the various cases 
that will be considered in the paper. The basic research question, relative to variable order 
operators, derives from the fractional differential equation 

c D?y(t) = f(t) (3.1.1) 

and the inferred integral equation 

c D~ q f(t ) = y(t ) (3.1.2) 

Since q in the fractional calculus can take on any real (or complex) value, the question is asked, 
“What is a desirable definition for the fractional integral when q is allowed to vary either with 
t or y ? ” More specifically, what is an appropriate definition for 

c D; q(t ' y] f{t) ? (3.1.3) 
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The analysis that follows is based on the Initialized Fractional Calculus (Lorenzo and Hartley, 
2000). The Riemann-Liouville basis for initialized fractional order integration is 

t 

c D; q f(t)s c d;‘ , f{t)+ii/(f-q,a,c,t)=-^j(t-r)‘ l ~'f(r)dr+y/(f-q,a.c,t). !>c, (3.1.4) 

<’ 

where c is typically taken as zero, and \j/{f -q.a.cj) is the initialization function. Only variation 
of <7 with t will be considered initially. Thus the problem becomes that of defining 0 D^ r) /(r), 
and determining the appropriate properties. A list of general properties, G7-G5, follows. 


3.2 Gl-Form and Bas is 

With the exception of section 4.5, consideration will be limited to fractional integrals whose 
forms are defined by 


o 



rfoM) 


f(j)dr 


(3.2.1) 


This eliminates from consideration, in this paper, such forms as 


-</(d 





r)" (, ' rl 

i , 

b 

rl 




fk)di 


(3.2.2) 


and many other possible forms which admittedly may be valuable. Some consideration of 
variable order based on the Grunwald form of the fractional differintegral may be found in 
Lorenzo and Hartley [1998]. 


3.3 G2 -Initialization of Variable Order Fractional Integrals 

We consider an initialization period a<r< 0, and require that /(/) = Ofor all t<a. For 

terminal initialization we require that, for ?>0, the initialized fractional integral starting at 
t - 0 be a continuation of the fractional integral starting at / = a, thus 

u D^'\f(t)= a D;^f(t) t> 0, (3.3.1) 

or in terms of the uninitialized fractional integral 

0 d;^f{t) +W {f,-q, a, 0.0= X“ /(,) /(') • (3-3.2) 


Therefore, we have 


' f ( f _ r y/(>.rM 

J r(q(t.r)) 


f(r)dr +y/(f, q{t.r\ a,0, r)= 


J r(</(r,r)) 




t> 0. 


(3.3.3) 


and 

>//■(/, q(t,z\ a,0. t)= J ~r f(r)dr , '>0- (3.3.4) 

a 

Thus it is seen that regardless of the argument of qit.x), the integrand of the initialization 
function will always be the same as that of the operator . 


3.4 Gi-Convolution Integral and Impulse Response Function 

It will be useful in the discussion that follows to 
draw on the ideas of convolution and impulse 
response as they relate to linear system theory. 

Briefly, we consider the block diagram of a linear 
system, figure 3.4-1, with input .v(r),and output 
y(r). The response of this system to a unit impulse Figure 3.4-1 System Block Diagram 
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function S (t ) is h(t) the impulse response function. The input, or forcing function, for this system 
may be written as an infinite sum of impulse functions, <5(/), 

.*(/)= J.v(r)5(/-r)</T . (3.4.1) 

Now for a linear system, superposition holds. Then, each of the input impulses in x(t) will excite 
a response in the output y(t) with a proportional strength and properly referenced in time to give 
the response 

y(t)= j x(r)h(t -z)dz . (3.4.2) 

A causal system only responds after an input is applied, thus /j(t)=0, fort <0. Furthermore, 
h(t -r)=0 for r >t , therefore we may express the response as the convolution integral 

t 

v(/)= j x(z)h{t-z)dz . (3.4.3) 

0 

Much more detailed discussion of these matters may be found in such texts as (DeRusso, Roy and 
Close). 

3.5 G4-Time Invariance 

From conventional linear system theory, with fixed and integer order, it is well known, (and 
important) that when the impulse response h(t, r)* hit -z) the system is said to be “time varying,” 
that is, the coefficients of the derivatives in the associated differential equations are functions of 
time (DeRusso). Fixed order operators (systems) are said to be time invariant if, when an input 
f(t) produces a response y(r), then the input /(r + cr) produces the response y{t+o). Because 
the result of variable order operation (differintegration) depends on two variables, the input /(/) 
and the order q{t ), this definition is inadequate for variable order operators. The concept will be 
generalized for the domain of variable order operators. 

A variable order operator will be said to be time invariant if an input / (?) and order profile 
q(t) produce a response y(t ) , then the input f(t + cr)and an order profile qit + o) will produce the 
response y(t+<y ) . 


3.6 G5-Linearit\ 

For the form above, equation (3.2.1), a// the operators are linear with respect to the input 
fit), that is. 


/ ( 'r m.t)) f(t)dr * b S : ^Xt)) (3 ' 61> 

It is noted that the operators are not in general linear with respect to the order q(t,r ), (i.e., q(t) as 
an input). Hence, we will not in general have the composition property, i.e.. 


.r)-l 




*r)-l 


,f(/-r) 9( ' 


■rU 


o d; u(,) 0 Dr i,] f(t) = 0 D; u ^/(t). ( 3 . 6 . 2 ) 

We will find that this is a peculiar aspect of variable order operators; it appears they can 
simultaneously be both linear and non-linear! 
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4. Toward Variable Order Operator Definitions: Variable Order Integration 


4.1 Kernels and Arguments 

We start by examining several Riemann-Liouville based variable order integrator definitions. 
For fixed fractional order integration we have for the uninitialized fractional integral 



c d t q f(t)= r /)j( f ^ f( T ) dT - 

C 

(4.1.1) 

When the order is allowed to vary with time, we write, taking c- 0 

f ('- T )‘ ?( '' rM f( T )dr 
01 n ) J r((/(/.r)) n > ' 

(4.1.2) 

or alternatively 

0 d , qU '/(/)= f -7 a/(f r)dr . 

Inq (t.r)} 

(4.1.3) 

where q* (t,r)=q(tj -r). 




In the discussion that follows we will look at the implications that result from various choices 
for the arguments of q(t.r). It appears reasonable to assume that the arguments of q in the 
exponent of (r-r)and in the gamma function, of equation (4.1.2) are the same, and that 
assumption will be made. The following cases will be considered; q(t,r)=q(t), q(t.r)= q(r ). and 
q(t,r)= q(t — r). One approach will be to consider a step change in q(t) from one constant value to 
another and to draw inferences relative to the character of the variable order integration operator 
based on the behavior of the kernel h(r.r), where 


h(t.r)~ 


r(<?M) 


(4.1.4) 


4.2 Case 1 q(t.z)=sq(t) 
4.2.1 Time invariance 
For this case we define 


o 


(t-r)« (,M 

r (q(r)) 


/(fV T. 


(4.2. 1.1) 


The associated impulse response here is 


h(t.z)-- 


(t-r )^ 1 

F(<?(')) 


(4.2. 1.2) 


Because the argument of q in the exponent of (z-r)is not a function of (t-z) it is seen that 
h(t, t)* h(t - r). Therefore, the associated dynamic system (the operator) is not time invariant. 

We now limit consideration to 0 < q < 1 , and write 




1 

r^Xz-r) 1 -^- 


(4.2. 1.3) 


To examine the behavior of this definition we consider a step in q(t) at t = , figure 4.2.1, and 

examine h(t y r)in the vicinity of the step. From this figure it can be seen that 0 d~ q ^f{t) under 
this definition immediately “forgets" that, from time t = 0 to time t = /, that the impulse response 
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(t — (t — rF 2 ” 1 

(kernel) was /i(?],t)= ^ — and switches with no “memory” of q t to h(t 2 ,z) = - ' 


rtffi) 

for t 2 > T m . . Specifically all “memory” of q x is lost after a change in q 


r(? 2 ) 



Figure 4.2. 1-1 Kernel behavior for q{t % x)-q{t) 
based definition 


4.2.2 Initialized Operator 

Now for starting time T z with t>T 2 > r w > 0 , and with f(t)= 0. Vr < 0 we have from the 
general property G2 


r; >/(>)= J ( f(r)dT+yf(f,- q(t\0.T 2 ,t). 

or 

t>T 2 . (4.2.2. l) 

r,D ' /(,)_ | rW» fiT)dr+ 1 rW>)) /(rVr ' 

t>T 2 . (4.2. 2.2) 

For the q(t) profile chosen this becomes 


r ' D ' ,0/ (' )= | ( n, s) /(t) * + | ( i4) f(t)dT • 

t>T 2 . (4.2.2.3) 


Clearly after the step change in q(t) all “memory” of q t is lost. That is, even the initialization 
function, which is the last integral in equation (4.2.2.3), does not “remember” q ] . In fact, all 
initialization is as if the operator (system) was always operating at the current value of q . 
Specifically, the properties of the operator shift instantaneously with q{t) and behave as though it 
has always been operating with the current value of q . 
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4.2.3 Physical Realization 
We consider q to be a general function of 
time, namely q = q(t ). We start by allowing 
q to change in a step-wise manner, more 
specifically, in a piece-wise constant manner 
that approximates q(t). The block diagram, 
figure 4.2.3- 1, may be used to represent this 
form of time varying fractional integration. 
Here, the switches, S ; - , associated with each 

constant q x integration, are normally open 
and are sequentially closed for duration T 
(initially finite) then opened again. In this 
realization, all of the fractional integrators 
start at time t = 0, run simultaneously and in 
isolation of the others. 

The following then applies 


Order 

Switch 

<h 

S, 

<h 

S; 


S3 





Figure 4.2.3- 1 Block diagram for physical 
realization of q(t,r) = ^(r) based definition 


Interv al Closed 
0 < f < T 
T < t < IT 

T <t <3T (4.2.3. l) 

( n — l)T <t< nT 


The order of fractional integration proceeds as q t ,q 2 ,q^ , ■■■ , that is, as a sequence of piecewise 



Figure 4. 2. 3-2 Approximation of q(r ) , / = r = 1 


constant-order fractional integrations. The q, approximate q(t) as indicated in figure 4.2. 3-2. The 
approach therefore is to fractionally integrate first /(/) at order q { , then at order q 2 , and so on. 
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The output of the process, y(t) , then is given by 


>’(') = 


Thus for any particular increment we have 

}'(') = o d ~ qk f(t) 



0 </< T 


d; qi f{t) 

T < t < IT 


d; q 'f{t) 

2T <t <3T . 

(4.23.2) 

d; q "f{t) 

(n-\)T<t<nT 



or 


>(')=/ 




Pk -i 


rfe*) 


■f(r)dr 


(k - \)T <t < kT, 


{k-\)T <t<kT. 


(4.2.33) 


(4.23.4) 


It can be seen from equation (4.23.2) that the final value of q is used at any point in time, t . 
Also, from the block diagram, figure 4.2.3- 1, it is apparent that during the time that any switch S, is 
closed, the output v(t) is as if q i had been in effect from the beginning. That is, 
q-q k , for ( k-l)T <t < kT , V& . We now take the limit of equation (4.23.4) as 7->0, then the 
stepped q k — > q(t), approach the original function q(t ), and on the left hand side of the equation, the 
output of the process y(t)->od?^ f(t ) . On the right hand side as T — » 0, n — > 00 such that for k = n, 
nT - 1 and we have as the definition for variable order integration for this type of process as 


(4.23.5) 



With the initialization included 


where 


0 D r -* w /(r)« 0 d; q(,) f(t)Mf -9(4 a.O.r), 

t> 0, 

ffr-rVW- 1 

</'(/• f(r)dr , 

a K 

r >0. 


(4.23.6) 

(4.23.7) 


The principal characteristic of this definition is its distorted “memory, ^ it always remembers the 
past as if the order of the operator q is at the current time, that is, it has no q memory. 

4.2.4 Operational Method q(t, r)=» q(t) 

Our interest here is to present a Laplace transform based procedure for dealing with 
variable order fractional operators based on the definition using q{t y r)— >q(t) i.e., those based on 
equation (4.23.5). All integrators (figure 4. 2. 3-1) operate continuously on the common input 
/(r). The result of the individual (fixed order) integrations we will call x(q i j). We observe 
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and 

L{x( qi j)}=r*F(s) (4. 2.4. 2) 

In the limiting process the number of integrations becomes an infinite continuum and the 
combined result may be viewed as a solution surface x(q(r\t) that is progressively sampled as 
q(t) increases. The solution then for the variable q(t)case is a path defined on this surface by q(t) 
which may be considered as a function on a new dimension /, . Because this definition has no 
memory of past values of q we consider two time bases; t as applies in equation (4.2.4. 1), and t x 
for the time varying order < 7 ( 7 ,). This allows us to describe the behavior of equation (4.2.4. 1) 


using an iterated Laplace transform. Specifically, we have 


(4.2.4. 3) 

(4. 2.4.4) 

In this formulation the iterated transform serves to order the procedure and after inverse 
transforming we will set r { -t . 

Example 

Consider the variable order fractional differential equation 

o d} a + b,) h*(t)=U{/). a >0, b>0, (4. 2. 4. 5) 

where it is desired to find the function h*(t) which when fractionally differentiated to the variable 
order a + bt yields the unit step function. Neglecting initialization, and remembering that 
composition does not apply we will instead solve a related integral equation. Using the dual time 
base, we write 


i ti (')]}= k L^'^'Cv)}. 

f-v, L {0 ^ (r ' ’/(')}= i-v, {u ?(,|) }f(4 0 < q(l)< m 


h(t)= 0 d; ia+b,,) U{t). (4. 2. 4. 6 ) 

Then transforming we have 

L Si L s M')}= L tl [l 5 ( 4 . 2 . 4 . 7 ) 

Using equation (4. 2.4.4) above 

L S] L, [h(t)}=-L Si L |ln, ' )1 ' ]- (4.2. 4.8) 

s s 

L s, L Mt)}=~ Si+ l^ • Re h (4.2. 4.9) 

Inverse transforming first with respect to s { yields 

L,{h(t)}=— e- b{]ns)t ' =l^L = s - (l+a+b ‘' ) , Re 5 ] >-Re(Mn(s)) (4.2.4.10) 
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Inverse transforming with respect to 5 , and replacing t x by / , yields 

a+bT { .a+bt 

h(t)=—. r — ^ r. l + < 

w r(i+a+^/,) r(i+<7+z?f) 

4.3 Case 2, q(t.r) = q(z) 

4.3.1 Time Invariance 

For this case the variable order fractional integral is defined as 

w-«(')/YA= f (?-r)‘ ?(r) ~ l f( \. 


l + a + bt > 0. 




The associated impulse response here is 

h(t,r)=- 


(/-r) ?(r> 


(4.3. 1.1) 


(4.3. 1.2) 


v ' r(q( r)) • V ' 

Because of the argument of q in both the exponent of (t-r) and in the gamma function it is seen 
that h(t,r)* h(r-r). Therefore, the associated dynamic system (the operator) is not time 
invariant. Limiting consideration to 0 < q < 1 , we write 


h(t,r)=- 


(4.3. 1.3) 


V ' r( ? (r)Xr-r) 1 -‘' (rr 

As in the previous case the behavior of this definition is studied by consideration of a step in q(t ), 
figure 4.3. 1-1, and examination of /i(r, r ) in the vicinity of the step at t = T m . 



i - integration variable 

Figure 4.3. 1-1 Kernel behavior for q(t,r) = q(z) based definition 


For q(t ) as shown in the figure, we have for t = t x < 

hAt, t)= — - — 0<T<t = t l . (4.3.1.4) 

r^X/.-r) 1 -^ 
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(4.3. 1.5) 


and for t = t 2 > T nv 


h 2 {t.t) = 


1 


i 




1-^2 


0 < t < T„, 2 

T**<*<t=h 


It can be seen that for t > T su . the impulse response function is segmented, i.e., it is 
discontinuous at as is q(t), This definition “remembers" the past q , namely q , in this case. It 
is also apparent that memory of ^ decreases (relatively) as t > 7^. increases. This behavior is a 
key difference when this definition is contrasted with the definition of Case 1. 


It is important to observe that under this definition 
both / and q are functions of r.i.e., /(r)and ^(r). 
When the operator is viewed as a dynamic system (block 
diagram of figure 4.3. 1 -2) it is not unreasonable to expect 
the “inputs" to have the same argument. It is interesting 
that we have here, and for all variable order operators, 
three terminal block diagram elements that explicitly 
contain neither an addition nor a multiplication. 


q(t) 



Figure 4.3. 1-2 Variable Order 
Operator block diagram 


4.3.2 Initialized Operator 

From property G2, and for a starting time T 2 > T m . > 0 and for /(/) = 0 Vr < 0 , we have 


T 2 D; q{,) f(t)=j 


Ht-ry {rU 


f(r)dr + y/(f,-q(t\o,T 2 j). t>T 2 (4.3.2. l) 


where y/(f,-q(r\0, 7\ . /)= j- 


T '( t-TY itU 
r (q(r)) 


,^2-1 


r(?(r)) 

f(r)dT . For the q(t) profile chosen this becomes 


rM 


■f(r)dr+ |- 


rh 2 ) 


-f(r)d t. t >T 2 . 


(4. 3. 2.2) 


The sum of the last two integrals here is the initialization function for T D, 


4.3.3 Physical Realization 

We again consider q to be a general function of time, namely q = q{t). Again, we start by 
allowing q to change in a step-wise manner, more specifically, in a piece-wise constant manner 
that approximates q(t). The block diagram, figure 4.3.3- 1, is used to represent this form of time 
varying fractional integration. The switches, S, , associated with each constant q t integration, are 
normally open and are sequentially closed for duration T (initially finite) then opened again. The 
switches in this case are on the input to the individual operators, and have the effect of parsing the 
dynamic input f(t). In this realization, all of the fractional integrators operate continuously from 
time r = Ountil timer that is, the end of the integration period. It is also noted that the 
fractional integrators continue to generate output even after the parsed input function f(t) 
becomes zero. The outputs are summed to form y(r). 
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Figure 4. 3. 3-1 Block diagram for physical 
realization of q(t,r)= <?(r)based definition 


The following then applies 


Order 

Switch 

Interval Closed 


Si 

0 <t< T 

2 

s 2 

T < r < IT 


s 3 

IT <t <3T 

<ln 

S„ 

(n~\)r < t < nT 


( 4 . 3 . 3 . 1 ) 


The order of fractional integration proceeds as that is, as a sequence of 

piecewise constant-order fractional integrations. The approximate q(t) as indicated in figure 
4.2. 3-2. The approach, therefore, is to fractionally integrate first f(t\u(t)-U(t-T)] at order q x , 
then f(t)[u(t-T)-U(t-2T)] at order <? 2 ,and soon. 

The output of the process, y(t ) , then is given by 

y(t) = o^‘ 9(r) l/(')((«0-«(f-7')))]+ 0 d~^ 2T * \f (t )((u(r -T)-u(t-2T )))]+ - + 

o d?> {nT) [f(t)((u(t - (n -l)T)-u(t - nr)))], (43.3.2) 

More briefly 

y(0=X o^r* (n ' ) KOW»-0-l)r)-4-»T)B. (4.33.3) 

/=] 
or 

y(0= X J~ f(^(iT)) ~ (« ~ W - «(t - iT)}]dx . (4.33.4) 
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Interchanging integration and summation gives 


(4 - V35) 

We now take the limit of this equation as n On the right hand side T — >0 such that 
t = nT. For clarity we introduce the dummy variable cr, and IT -+o , and 

f(t)^f(iT)= f(a), q(iT )— » qip) and the quantity ^ r — — — > 5{r-o) the 

impulse function and we have 


unit 


*> II 

0 0 




-/(cr)d(r -o)do dr. 


(4.3. 3. 6) 


Through the limiting process the left hand side becomes the desired variable order operator and 
we have the definition for this case as 



(r-r )^ 1 

r(<?(r)) 


f(t)d r. 


(4.33.7) 


4.3.4 Operational Method 

Here we seek an operational method for solving equations (systems) with variable order 
fractional integrals based on the definition 


»‘ , ’ < '’ /( ' ) *| i Trar /(r) '' r (4 -- vu) 

It is useful to refer the block diagram, figure 4.3.3- 1 . Conceptually we will take the finite array of 
fixed order integrators of order q , J = l,2,...n to the limit to become an infinite continuum of 
integrators. That is, an infinite array of fixed order fractional integrators of order q, - 4(A) where 
A is allowed to take on all real values 0< A <«!<<*>. We now convert /(f) into a A distributed 
input, thus we define 

.v(A,/)=/(A)5(f-A), (4.3.4.2) 

as the input directly into the integrator (after the switches), and observe 


x(A, f) = 



A = r 
A = 0 


(4.3.43) 


We further note that 


f(t)= J.v(A,f)5(A-fVA . 

^0 

The output of the individual integrations, ©(A,f), is 

0(A,f)= o ^ <A) .v(A,f), 

and the integration over A of these outputs will be 

Q d;^f(t)= J e(A,f) dX = j 0 d^ x) x(x,t)dX , 

^0 ^0 


(43.4.4) 

(4. 3.4.5) 
(43.4.6) 
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(43.4.7) 


4* 

0 </,- ?(,) /(')= J 0 d;^ ) f{X)8{t-X)dX. 

Xq 

The Laplace transform with respect to t , with Laplace variable s , then is 

L,{o^V(7)}= W ]od; qW f{t)8{t-X)d\ J = ~je~" J 0 d^ {x) f(X)S(t-X)dXdt. (43.4.8) 

Ao J 0 Aq 

Under appropriate conditions the order of integration is interchanged giving 

^{o^ (,) /(')}= jf(X)je- s ' 0 d;^ ) 8{,-X)dtdA, (43.4.9) 

0 

^{oV (/> /(0}= ^mr^e-^dk. (43.4.10) 

Ao 

At this point there are two possible paths, the first is to recognize that A is a secondary time 
base that we will call t { and write the time varying Laplace transform 

(4.3.4.11) 

In application of equation (4.3.4. 1 1 ) the order of solution is; inverse transform with respect to s , 
integrate over A , and replace /, by t . 

The second approach, is to take the Laplace transform with respect to A of equation (4.3.4.10), 
taking A () = 0 , this gives 



L-s, L {o ^ ?(,) /(')]}= 4- L s> {f (A)^< V* }, (4.3.4. 1 2) 

(4.3.4.13) 

In application of equation (4.3.4.13) the order of solution is; inverse transform with respect to s , 
integrate over A , inverse transform with respect to s, and replace /, by t . 

Example 

We now consider the same variable order fractional integral equation (4.2.4. 6) of the previous 
case, however now with the definition based on q{t,x)=q{t). Thus we seek to find the function 
/i(/), i.e., which is the variable (a + />f) order fractional integral of the unit step function. We 
neglect initialization, and using the dual time base, we write the integral equation. 



Then applying equation (4.3.4. 1 1) we have 


L s (f)}= L s {/?(?)}= jtf(A)s~ (a+a) e _/U </A 

0 

e a 4 


(4.3.4.14) 

(4.3.4.15) 

(4.3.4.16) 
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| _ e U+b\n(s)) 

.V + frln(5) 


(4.3.4.17) 


This result should be compared with that of the previous case equation (4.2.4.10). the forms are 
found to be substantially different. Then the solution is given by 


h(t)= 



ln (d) 

s"(.v + Mn(.9)) 


} h f 


(4.3.4.18) 


4.4 Case 3. q(t,r)~ q(t -v) 

4.4.1 Time Invariance 

For this case we consider the definition 


jVJ,-,)) /<rVr - 

( 4 . 4 . 1 . 1 ) 

The associated impulse response here is 

r(?(f-r» 

( 4 . 4 . 1 . 2 ) 


Because the argument of q in the exponent of (;-r)and in the gamma function is of the form 
(r-r)it is seen that /?(r,r) = /i(r-r). Therefore, the associated dynamic system (the operator) for 
this case is time invariant , as defined in property G4 earlier. 


We now limit consideration to 0 < q < 1 , and we write 

1 


h(i, t)=- 


-. (4.4. 1.3) 

r(g(f-r)Xr-r)'- ?, '- r) 

We again examine the behavior of this definition by considering a step in < 7 ( 1 ), figure 4.4. 1-1, 
and examining /j(r, r ) in the vicinity of the step at / = T m . For q(t) as shown in the figure, we have 
for i = t t < T w 


h\ ('•*■)= 


r(<?,Xb-t) 


l-q ] 


0 < T < t = t x . 


(4.4. 1.4) 


and for t = t^ > f 


/? 2 (r,r)= i 


1 


1 

r fairs 


0 < r < t - T n 


T„. < r < r = r. 


(4.4. 1.5) 


As in Case 2, it is seen that the impulse response function is segmented, i.e., it is discontinuous 
at This definition very strongly “remembers” the past q, namely q^ in this case. This 

effect is so strong thal if we consider the integral at time t = t 2 where t 2 =7^. +e , fore small, we 
have the situation where even though q has switched from < 7 , to q 2 the variable order integral is 
almost completely evaluated at q { ! This is a key difference when this definition is contrasted 
with the definitions of the previous cases. 
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x - integration variable 


Figure 4.4. 1-1 Kernel behavior for q(t,r)= q{t-r) based definition 


4.4.2 Physical Realization 

We again consider q to be a general 
function of time, namely q = <?(/). Again, we 
start by allowing q to change in a step-wise 


manner, more specifically, in a piece-wise 
constant manner that approximates q{t). The 


block diagram, figure 4.4.2- 1, is used to 
represent this form of time varying fractional 
integration. The switches, S i , are now each 


associated with a constant q n _ { integration. 


It is seen that the sequence of integration 
now runs backward relative to the switching 
of the input function f(t). This necessitates 
knowing the entire q(t) profile in advance. 


The switches are normally open and are 
sequentially closed for duration T (initially 
finite) then opened again. 



Figure 4.4.2- 1 Block diagram for physical 
realization of <?(f,r)= ^(r-r)based definition 


The switches in this case are on the input to the individual operators, and have the effect of 
parsing the dynamic input f(r). In this realization, all of the fractional integrators operate 
continuously from time t = 0 until time t-t, that is, the end of the integration period. It is also 
noted that the fractional integrators continue to generate output even after the parsed input 
function /(/) becomes zero. The outputs are summed to form v(r). 
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The following then applies 


Order 

Switch 

Interval Closed 

Qn 

S, 

0 <t< T 

tfn - 1 

s 2 

T <t <2T 


S, 

2 T <t< 3 r 

q : 

S„., 

(/2 - 2)T <r<(»- 


S n 

(/? - i)r < r < /7r 


The order of fractional integration proceeds as q„ ,q „_ 2 , that is, as a reversed sequence 
of piecewise constant-order fractional integrations. The q : approximate q(r). The approach 
therefore is to fractionally integrate first /(/)[f/(r)-£/(/-7’)] at order q„, then 
f(t)[u(t-T)-U(t-2T)] at order q „_ { , and so on. 


The output of the process, y(r ) . then is given by 

y(r)= 0 ^ (,,r V(0(W-^-7')))]+ 0 ^«''-‘> r >|/(,)(( M (r-r)- M (r-27')))]+ - 

More briefly 






i-0 


or 


><')= u 




“f{ r(q((n-i)r)) 

Interchanging integration and summation gives 


[/'(r){w(r - iT )- u(z - (/ + 1 )T )}]</z 


(4.4. 2.2) 


(4.4.13) 


(4. 4. 2.4) 


*>-js 




r(q((n-i)r)) 


/( T) 


{t/(r-iT)-n(r-(/ + l)r)} 


T\d x . 


(4.4. 2.5) 


We now take the limit of this equation as n — > °o . On the right hand side T — > 0 such that 
t = nT , for clarity we introduce the dummy variable a , and iT ->o , T — » da , and 

f(z)—>f(iT)= /(<r), qijT ) q(o) and the quantity M r ~ lT U ( T ~ ( f + -> <5( T - cr ) the unit 

impulse function and we have 




(4. 4.2. 6) 
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Through the limiting process the left hand side becomes the desired variable order operator and 
we have the definition for this case as 


o 




f(r)dT. 


(4.4. 2. 7) 


It is apparent from this realization that the (r -r)argument of q effectively causes <7 to run 
backwards in r ; that is, the entire q(t) profile must be known apriori. Thus, for this definition, 
q(t) appears to behave as an integral part of the kernel not as an “input” to the operator. 


4.4.3 Initialized Operator 

Directly applying the G2 property, for the starting time t - 0 > a , and for f(t)= 0, Vr < a 
the initialized operator is given by 


. = 0 d^fit) + y/(f-q(t ), a,0, t ), 


t > 0, (4.4.3. l) 


or 


r>0 - < 4A3J > 

0 a 

Then for the test condition associated with figure 4.4. 1-1 with the starting point of the integration 
at T 2 > T > 0 , and for /(;)= 0, Vr < 0 , we have for r > T 2 


-r)-l 




■■ D '"' ,f{,) * I iCr /<r) JT + ’ f *?$r*»**j /(rVr • (4A3 ' 3> 


where the sum of the last two integrals is the initialization function. 


4.4.4 Laplace Transform 

The derivation of the Laplace transform of the variable order integral follows that for the fixed 
order case, since the convolution theorem can be applied. Then, considering the uninitialized 
case, i.e., Q d;‘ ,(r) /(r)= 




pr, <7(r)>0,r >0, (4.4.4. l) 


then the Laplace convolution theorem is given by 


L{u(t)*v(t)}=U{s)v(s)=l] 






(4. 4.4. 2) 


fHV 

Now taking u{t)= f{t) and v(t)= . . . . = h(t ) the convolution theorem yields 

n?w) 

4 0 D; qi,) f(i)}= u(s)v{s)= l{/(/)}l|-^L|. 


(4.4. 4.3) 


It is important to observe that under this definition q is fully associated with the kernel impulse 
function i.e., q-q(t- r), that is, the argument of q is the same as the (t-r) term that 
it exponentiates. When the operator is viewed as a dynamic system, block diagram of 
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figure 4.3. 1-2, it is not unreasonable to expect the “inputs” to have the same argument. This is the 
case when q(t.r) = q(z) but not so here. 


4.5 Alternate Approach Variable. Order Operator Based on Laplace Tnmsfonn 

The following approach uses an extension of the Laplace transform for fixed order operators to 
infer a possible definition for the variable order operator. Neglecting initialization (or initial 
condition terms), and assuming functions of exponential order, we note the following 


J- J/(LVr„-dt =5-' , F(.v) 



= s~'F(s) 


L {il r ( ,)j = .,= f ( s) 


(4.5.1) 


Thus, in general we have 

4 >^/(')}= s?/r (*)’ v <7- ( 4 - 5 - 2) 

Now since Z.{ct/(/)}= c/s , equation (4.5.2) may also be written as 

V*. (4.5.3) 

In this form <7 is of course considered as a constant. However, it is now a simple matter to 
generalize q by considering it to be a function of time. Therefore, a Laplace transform based 
definition for non-initialized variable order integration and differentiation is 

(4.5.4) 

If it is assumed that the order requires no initialization then the definition for the initialized 
operator is given by 

t> 0, (4.5.5) 

where we require 0 df'^f(t) and q(t) to be of exponential order and piecewise regular. While this 
unorthodox definition lacks a clear time domain meaning and realization, it is operationally 
compelling and may find important application. 
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5. Characterization of Variable Order Fractional Integrals 


5.1 Generalizing the a- Argument 

Through the study of the previous cases we have seen the behavior of three time-based 
possible definitions for the variable order fractional integral operator. We have seen that the 
operator based on q(t,z)=q(t) has no memory of past q , the operator based on <y(/,r ) = ^(r)has 
weak memory of past q , and the operator based on q(t,z)~ q(i -z) has very strong memory of 
past q . It is clear that the response of the operator to changes in q will be inversely related to its 
memory of past q. These qualitative results are summarized in the table below, where a 
generalized operator is defined based on 

q{t.z)=q{at+bz). (S.l.l) 

that is 


°d' q{ ' )f ^~^r(qlr + bz)) fij)d 


[a r+/?r)-l 


(5.1.2) 


Case 

a 

b 

Relative q Memory 

Response to q Change 

i 

i 

0 

none 

immediate 

2 

0 

1 

weak 

intermediate 

3 

1 

-1 

strong 

very slow 


This summary suggests that the behavior of a variable order fractional integral operator might 
be '‘tailored” to a particular application by suitable selection of the constants a and b in the 
generalized form, equation (5.1.2). The rest of this section will create some tools that will be used 
to evaluate such definitions and to characterize their behavior. The first task is to determine the 
useable range for selection of a and b . 

For the defining uninitialized operator, equation (5.1.2) above, the variable of integration r , 
ranges as 0 < r < t . Therefore, we require that the argument of q also fall in this range, i.e., 

0 <at+bx<t. (5.1.3) 

If the lower limit is violated q is not defined, violating the upper limit the operator is non-causal. 


Examining the boundaries, from t-x we 
have 0<(a + 6)r<f or 0<a + b<\. From the 
condition r=0, we have, 0<af + 0</ or 
0 < a < 1 which together with the previous 
condition yields -1 <b<\, constraints that are 
preempted. The acceptable range of a and b 
is shown in figure 5.1-1. Also shown on the 
figure are the three cases analyzed previously. 
It is noted that for certain non-real time 
applications, such as digital filtering, it may be 
desirable to relax the upper constraint, 
i.e., causality, and allow the use the future 
values of q . This has the effect of removing the 
upper and right constraints in figure 5.1-1. 


b 



Figure 5.1-1 Allowable range 
of q(at+bx ) argument 
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5.2 Memories 

Two types of memory are considered. The first type, m { , is the memory associated with the 
integral operator itself. This is the so-called fading memory effect of the fractional integral when 
compared to the integer order, <7 = 1, integral. A second type of memory, m 2 , is that associated 
with the change of order as we have indicated above. In the material that follows, the focus will 
be on variable integrators of order q(t) bounded as 0< q{t)^ 1 . 

5.2.1 m { memory 1 

This type of memory is the well understood fading memory of the fractional operation and 
will be discussed only briefly. We first consider the fixed order fractional integral 


C7(')= J 


H± 


rT‘ 

r(<?) 


f(*)d 


T. 


(5.2.1. 1 ) 


In this integral the past is described mathematically by (t-r)in the kernel. When r is small or 
zero we are at or near the present, i.e., / -r — > t , conversely when z is large or approaches t in 
magnitude we are in the distant past, i.e., t-x ->0. Thus we integrate, fractionally, from the 
present, r = 0 or r — r = f, to some point in the past, r = /or r-r= 0. The form of the kernel, 
through q in the exponent and in the gamma function, determines how the past information is 
weighted. The response of the integral for a unit impulse function at time / = 0 is 


7<r)-J 






S(r)dz =^—r-^- = h(q,t,v = 0 ) . 


(5.2. 1.2) 


It?) V ' It?) 

u 

which is the operator impulse response function. The effect of the fading memory may be 
observed by comparing the response for fractional q , say after t = 1 , to the integer order case 

<7 = 1, which remains constant at a value of unity. 


The impulse response function may, 
therefore, be used as an indicator of 
m x type memory. Thus, a primitive 
measure of the memory of the fixed 
order integral operator is the shifted 
response 


m 


< 0 - 


= h{qj\ 


' (fir 

. r w Jal 
0<<7<1. ?>1. 


f>i 


(5.2. 1.3) 


This measure is shown for the fixed 
order fractional integral in figure 5.2.1- 
1. When the order is variable with r, 
(time) this measure shows a variation 
with q(t). The following is an 
alternative measure for the m { type of 
memory 



Figure 5.2. 1-1 m x - Memory measure 
vs. time for fixed order fractional integrals 
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Here we observe, for f(t) = u(t) a unit step at t = 0 , that the ratio 



1 

r(<?(? . t )) 


7/ 


dr , 


0< <?(r,r)< 1, 


(5.2. 1.4) 


0 

is essentially the same as expression (5. 2. 1.2) above when q is of fixed order. This measure 
also shows variation when q is a function of time. These modified measures will show the 
fading memory effect with variable q . Better measures of m { are certainly possible, but are 
outside the scope of this paper. 


5.2.2 nu memory 

A more important memory type for the study of the “tailored” variable order fractional integral 
is the m 2 type that is associated with memory of the order history. While the m, memory is meant 
to be a measure of the retentiveness of the entire fractional integral, here we seek a measure of the 
retentiveness of the order within the fractional integral after a change in order. Thus we consider 
a step change in order from q x to q 2 occurring at time t-T as shown in figures 4.2-1, 4.3. 1-1, 
and 4.4. 1-1. Setting /(/)=! eliminates the effects of changes in the input function. We define m 2 
memory retentiveness as 




01 - 1 




rfei) 


dr 


area under q x portion of 


[M^ 
J r(<r(t,T 


M 


0<q(r 7 r)< 1, (5.2.2. l) 


area under h(t, r) 

„ r(,(,.T)) ' dr 

where x^and r 2 are the lower and upper bounds respectively associated with area under the q x 
portion of /?(r, r). Thus m 2 measures the strength or retentiveness of the memory of q x at time t 
when t > 7\ . 


This measure may be readily applied to the variable order definitions of Cases 1 to 3. It is 
readily seen that for Case 1, q(t,r)= q(t), (refer to figure 4.2. 1-1) 

« 2 (')= t ( 5 - 2 - 2 - 2 ) 


(o t>T m . 

Then for Case 2 above, ^(r,r)= q(r), (refer to figure 4.3. 1-1) 


''{t-T) 


l^l-l 


m 2 (f)= — 


r(?,) 


-dr 


(t-T) 


Q \ - 1 


r(?,) 


-dr + 


f (t-rr 


J 






-dr 


I > , 


(5.2.2.3) 


and for Case 3, q(t,t)= q(t -r), (refer to figure 4.4. 1 - 1 ) 
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(5.2.2.4) 


« 2 (0=- 


T m , 


t-T n . 


(t-t) 




rfo.) 


-dr 


i*7i “I 


f iziL 

r(<7.) 


/-7v 


f > L. . 




-</r + 


;-r„. 


fM 

i rfe 2 ) 


-rfr 


The behavior of the w 2 -measure for these three cases is shown graphically in figure 5. 2. 2-1, and 
thus the results discussed subjectively above, for these cases, have been quantified. 



Figure 5. 2. 2-1 m 2 (t)v s. / -Time for =1, 
q(t,T)=q(t), = q{ri = q(t - r) 



Figure 5.2. 2-2a /?/->( 3) memory behavior for 
‘Tailored" variable order fractional integrals 


rri^p* 



Figure 5.2.2-2b Contours of constant /w 2 (3) memory 
for “tailored" variable order fractional integrals 


To examine the m 2 behavior for the “tailored" integrals, equation (5.1.2), we select a 
particular value of t , say t =3 , in the m 2 -measure. This result m 2 ( 3) versus a and b is shown in 
figure 5.2.2-2a. This figure shows the variation in order retention (m 2 memory) as the a and b 
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parameters of the order q(a t + b r) of the variable order fractional integral vary over the allowable 
range. Figure 5.2.2-2b presents the same data as figure 5.2.2-2a but now as contours indicating 
loci of constant values of m 2 ( 3). It is important to note that the contours are straight lines 
emanating from the point (l//,o)=(l/3,0). From such data one may select the desired properties 
for the variable order integrals (and derivatives) for a specific application. 

53 Dynamic Behavior: Response to q Changes 

The responsiveness of the variable order fractional integral to changes in the order q is an 
important consideration in the application of these operators. In general, the longer the memory of 
past q , i.e., the greater its m 2 retentiveness, the slower will be the response of a particular 
definition to a change in q . This effect as well as the responsiveness of the variable order 
fractional integral to changes in the input function /(/) will be illustrated in the application 
section. 

6. Variable Order Differentiation 

The above definitions for variable order fractional integration may be formally extended to 
variable order fractional differentiation in a manner similar to that done for constant q . For the 

uninitialized case we define 

0 ^ f) /(/) = 0<0 d; u(,) f{t) t> 0, q{t)> 0, V/, (6.1) 

where mis an integer such that, q{t)= m-u(t\ and m><?(r)>0, V/. Similarly for the initialized 
case we formally define 


0 D;<'V('K*>r oDr {,) f(t\ t> 0 , q(t)> 0 , ( 6 . 2 ) 

where q(t)= m -u(t) l and m is taken as the least integer greater than the maximum value of q(t). 

Composition is inferred in the above definitions. However, the determination of a general 
composition law for the variable order fractional integral may subject these definitions to later 
reconsideration. Lack of a composition law increases the difficulty of analytical studies using the 
definitions for both the variable order fractional integral and derivative, however, numerical 
approaches are viable. Another complication is the fact that it may be desirable to allow <?(/)to 
range over both positive and negative values. This places an analytical “seam” at q - 0, requiring 
the use of different definitions at each crossing, which may cause difficulties when using an 
approach based on the Riemann-Liouville definition. Thus an approach based on the Grunwald 
definition may be required. A consideration of the variable order differintegral based on the 
Grunwald definition is presented in Lorenzo and Hartley, 1998. 

The character of the variable order fractional differentiator will depend on the character of the 
basis definition used for variable order fractional integration. That is, the memory and response 
behaviors of the variable order fractional differentiator will mimic those of the basis variable 
order fractional integrator definition. Thus for example, using the Case 1 definition q{t,x)~^ q{t), 
the fractional variable order derivative will have no memory of past order q and will behave in a 
manner similar to the integer order derivative. 
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7. Application of V ariable Order Operators 

The motivation section of this paper presents a number of potential application areas for the 
variable order operators. The introduction and consideration of the “tailored” variable order 
operators also allows application to a variety of adaptive filtering, detection, and modeling 
applications. The scope of this paper does not allow detailed study of these. This section will 
examine some simple applications to show some of the behavior of these definitions. To study the 
response of the variable order fractional integral to a change in q{t) we start by holding /(/) = 1 . 
Figure 7-1 shows numerical results for the three basic cases studied earlier to a square wave 
“input” of q{t). The nu_ memory retention of past q is readily observed by comparison of the three 

responses. For example, for the q(t,z) — » <?(/) transient the switching between the constant 
<7 = 0.75 and q = 0.25 responses is obvious. A more detailed view of the response to a step 
change in q(t) is shown in Figure 7-2. Here attention is drawn to the response of the integral after 
time / = 0.1 , the effects of the q argument in the kernel become very apparent for the three cases. 

We now consider a variable order integration with ^(r)= 0.5+ 0.35 sin(lO^t); then we have 

j -(0.5+0.35 sin (1 Off r ))y A \ f ( / r , (7. l) 

0</ ' J r(q(t,r)) 

where again the three definitions of ^(t.r)are applied. Figure 7-3 shows the results of these 
variable order fractional integrations. It is observed that the Case 1 definition q(/,r)— » q(r) result 
starts out 180 degrees out of phase with the “input” of q(t) but then shifts to an in-phase response, 
the Case 2 definition q{t.z)- > q(r) result remains 180 degrees out of phase with the “input” of 
q(t) for the entire transient, while for the Case 3 definition q(t,z)^> qit - r) the sinusoidal effect 
of q(t ) is barely perceptible in the result. 

The effect of the integrand is now considered. Figure 7-4 shows the results for /(;)= sin(6fff) 
while maintaining q(r)= 0.5 + 0.35 sin(lO;r r). The Case 1 definition q(t,r) — > q(t) and the Case 2 
definition q(t.z)-> q(z) results appear to be quite similar and have apparent frequency content at 
the sum frequency sin(l6^t). Close examination of the difference of these two results however 
shows content at the product frequency (not shown). For the Case 3 definition q(t,r)—>q(t-T), 
the result follows the integrand with little, if any, phase shift. 

That the integrand can interact with the order “input” is demonstrated in figure 7-5. Here again 
the Case 1 definition q(t,r)->q(t) and the Case 2 definition q(t,z)-* q(r) results appear to be 
quite similar, both having significant harmonic content. The startling feature of this figure 
however is the amplitude growth observed for the Case 3 definition q(t, r) — > q(t — r ) result. 

The matching of the frequencies, apparently, allows energy to accumulate in the resultant 
output of the integration. This is a resonant-like condition. It is also observed in figure 7-5 that 
the q(i,x)=q{r) result is slowly ramping up, also an indication of energy accumulation. Thus, the 
physical inference may be drawn that: varying the order requires (or yields) external energy. 
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t - Time 

Figure 7-2 Variable order fractional integral responses 
for arguments q(t,T) = q(t),= q(r), = q(t-r), q(t) input is a step at t = 0.1 
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Figure 7-5 Variable order fractional integral responses for arguments 
q(t,r) = <?(/),= q{r) s - q{t -x) , q{t)= 0.5 + 0.35 sin(l0;r f), /(/)= sin(l0/r /) 


8. Distributed Order Operators 

8.1 A pproaches to Distributed Order Operators 

As in the case of time variable operators, because the order of the fractional integral and 
derivative may take on any real value (or imaginary) we are allowed to consider variations of the 
order with some variables other than time, say for example, space. This section explores the 
distributed order operator. Two approaches are considered: a direct approach that does not assign 
a new variable for the variation in order and an independent variable approach that considers q to 
be a function of some independent (perhaps spatial) variable. 

8.2 Direct Approach 

The approach to be discussed here is a general mathematical method. However, to aid 
understanding we will develop the concept around a reference mechanical system. The general 
form for an element in such a system is 

*0 d?y(t)=f(t), (8.2.1) 

where, q = 2 for an inertial mass, q = 1 for a dashpot, q = 0 for a spring, y(t) is the displacement 
(response), and /(/) is the time varying force on the component. Other possible elements may 
include 0< q < 1 for various viscoelastic materials, and 1 < q < 2 for viscoinertial elements. Thus 
an assembly of such elements in parallel would be represented by 

X k ‘ (o d ? ><') + V'/ )= 2 fi (')’ (8-2-2) 

i / 

where i//, =\ff{y,q jy a,0,t) is the associated initialization for each element. 
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Some materials have rheological properties that depend on temperature, electrostatic field 
strength [Shiga] or magnetic field strength [Davis]. Bagley (1991) discusses materials that display 
complex thermorheological behavior. This means that the order of the viscoelastic element, q n , 
depends upon the temperature of the material. If a sample of this material is subject to a 
temperature distribution, a corresponding order-distribution will exist throughout the material, as 
each individual element will have its own order. In the limit as the elements approach differential 
size this leads to the concept of the order-distribution. While this concept has been introduced 
using a mechanical example, it is conceptually simple to jump from the summation in equation 
(8.2.2) to the integral in equation (8.2.3) below. Assuming elements with orders that vary from 
zero to two, the general mechanical system of equation (8.2.2), is replaced by an integral over the 
system order, 


J*(?) 0 d? y(t)dq + J k(q)y/(y\ q, ,0, t )dq = /(f), (8.2.3) 

o o 

where the first integral is named the cumulative order distribution over the range 0< q < 2 . With 
k(q) being called the order-weighting distribution or the order strength distribution. In general we 
will use the notation 

<1: 

„,. 9: J k(q) e d?y(f)dq , (8.2.4) 

< 7 1 

for the cumulative order distribution, where the q in the argument of £1 is the integration 
variable. For the special limits this will be written as 

O?, (k(q\ v(t\ q) = j k{q) c d?y(,)dq . (8.2.5) 

0 

For convergence of the integral (8.2.3) we shall require that the integrand be bounded and that 
A-(<?)be non-zero only over a finite range. The second integral of equation (8.2.3) is the 
cumulative order distribution initialization, and the forces have been combined into a single term. 
The general notation for the cumulative order distribution initialization is 

(*(?). v('W)= jk(q)y(y,q,a,c,t)dq, (8.2.6) 

<7i 

The Laplace transform of equation (8.2.3), assuming the integral converges and is of 
exponential order, is given by 


/( 4 )s v dq 


yu>+ 


i j k(q)L(y{y, q.a.Oj ))dq 


F{s). 


(8.2.7) 


Thus we now have a very general formulation for representing dynamic systems. Indeed the sum 
of all orders of constant coefficient linear derivatives can be expressed as a single term! 


To demonstrate that familiar equations can be written in this form, the common second order 
mechanical system may be readily expressed using equation (8.2.7). Setting the initialization 
integral to zero (or including it with the force), gives 


j[m8(q — 2) + b8(q-\) + dq 


\y (s) = f(s ) 


( 8 . 2 . 8 ) 
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or 


(ms 2 + bs + k)Y(s) = F(s). 


(8.2.9) 


The range of system order may be extended in both directions, then the general system 
representation becomes 


n l 
Y(s) + 


J *(<?)•*’ dq 

which has the time domain representation 
JV'(<7) o Dfy(t) dq = f(t). 


jk(q)L(y/(y.q,a.O,t))dq 


= F(s) 


( 8 . 2 . 10 ) 


— oo < q < < 


( 8 . 2 . 11 ) 


and k(q ) is constrained to be non-zero only over a finite range to assure convergence of equation 
(8.2.1 1). The forcing term on the right hand side of this equation may also be a cumulative order 
distribution, giving the more general form. 


f oo N 


/ oo > 


( °° 

dq 

Y(s) = 

J g{p)x p dp 

k _ °° J 

F(s) — 

1 


J k(q)L(\f/(\\q, a,0, t))dq 


( 8 . 2 . 12 ) 


This form contains all non-partial constant coefficient differential equations (both integer 
order and fractional) when the lower integral limits are set to zero! With the limits as written, it 
also includes all possible integral terms. Further extension to include time varying coefficients 
may be obtained by replacing k(q) with k(q,t) in equation (8.2.11) and transforming 
appropriately. The time domain solution is given by 


*)=i- 


°jg(p)s p dp 

. _ / _l , 

< -- r \s ) 

jk(q\x q dq 

* — s 


J k(q) L(y/(y,q,a,0,t))dq 
J k(q)s q dq 


(8.2.13) 


The analysis section considers some specific systems that can be represented by equation (8.2.10) 

8.3 Independent Variable Approach 
8.3.1 Definition 

This approach to distributed order operators considers q as a function of some independent 
variable x , (here considered to be spatial). Thus we have q -q{x) and k = k(x) and we consider 
operators of the form 

o df x) f{t). (8.3.1.1) 

For this differintegration q(x) is considered to be constant in the fractional operator (relative to 
the t (time) based fractional integration). Operators of this form will be explored in more detail 
in section 10 of this paper. The cumulative spatial order distribution then is defined as 

B 

A . B nff\k(q\f(t\ X )^jk(x) c d q{x) f(t)dx. (8.3.1.2) 
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Normally we shall take the limits of the .v integration as 0 and °° , and assure convergence of the 
integral by limiting k(x) to be non-zero only over a finite range of .v . Taking c=0,we will write 
for the cumulative spatial order distribution 

©o 

nfM‘l\f(t\x)=\k{x) 0 df' ) f{t)dx, (8.3.1.3) 

0 

In this definition ^(.v) is considered to be constant in the fractional operator (relative to the time 
based fractional differintegration). The general notation for the cumulative spatial order 
distribution initialization is 


B 

A.k'vfc'ikicil fitly) = J k(x)y/(f, q , a, c. /, x)dx , (8.3. 1 .4) 

A 

where the limits of the integral may be modified to suit the problem to be solved. From an 
engineering or scientific view .v might be the distance along some thermorheologically complex 
material that is experiencing a temperature gradient. Thus at each point .v the material responds 
with a different order q(x) and spatial strength spectrum k(. v). The application (figure 9-1 ) of the 
next section shows such a situation. Extension to multi-dimensional cumulative ( spatial ) order 
distributions is obvious 

.. V2 

J- * J — Jf„) V(r) , (8.3. 1.5) 

0 0 0 

where again the limits may modified to suit the problem to be solved, and k is constrained to 
assure convergence of the multiple integrations. 

8.3.2 Laplace Transform 

The Laplace transform of the cumulative order distribution (8.3. 1.3), assuming that it is of 
exponential order, is 

L^lf [x \k{q\ f(t\ ,)}= J J k(x) o d? {x) f(t) dx 

[o 

= J e- J k( X ) 0 df' ] f(t) dx dt = nf x) (k(q\ f(s\ .x). (8.3.2. 1) 

o o 

Under conditions allowing the interchanging the order of integration, we have 

= J o df x] f{t) dt dx , (S.3.2.2) 

0 0 
or 

J k(x) Q d^f(t) dx = J A-(,)l { 0 df x) f(l)}dx . (8. 3.2.3) 

0 J 0 


j*(.v) o df x) fit) 


dx 
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Since q(x) is considered to be constant in the fractional operator, this becomes 

| oo I oo 

f (t\ *)}= Lj J k(x) 0 df'\f(t) dx l = J k(x)s q(x) dx F(s) . (8. 3.2.4) 

In application it will be often desired to also transform in the x dimension, thus using s, as the 
Laplace parameter with respect to j: we have, limiting the integration range to v 


L S[ 


L S jj 40 odf x \f{t) rfvj | = L Si j J k(x)s qU) dx F(v)| 
= — L 5| t(-v)^ (v) }F(0. 


(8. 3. 2.5) 


(8. 3.2. 6) 


It is interesting to compare this result to the variable order results equations (4.3.4.12) and 
(4.2.4.4). 


9. Analysis with Distributed Order Operators 


It is useful to reconsider the transform of the cumulative order distribution, equation (8.2.5), 
by rewriting the exponential in s . That is 

j DO j OO OO 

44,(4<0 y(0<?)}= 4 /*(*) <d?y(t)d<i = J 4<?>^0V<7 = J4<7> 9ln(l V(0, ( 9 . 1 ) 

lo Jo 0 

This cumulative order distribution transform is effectively a Laplace transform of the function 
k{q) with the new Laplace variable being r = — ln(jr) . As long as the order-distribution, k(q) is 
such that the integral converges, the integral is easy to calculate by treating it as a Laplace 
transform with Laplace variable (-ln(s)). 


Consider the direct (approach) order distribution k(q) = U(q)-U(q- 2), that is, a uniform 
distribution of magnitude 1, between q- 0, and q = 2. The Laplace transform of k(q) from the 
tables is 


r r 

Then replacing r by -ln(.v) gives 


oo 2 

^k(q)e qUs) dq =^e Ms) dq = 
0 0 

for the transform of the distribution. 


i 2 l n(s) i 2 

1 -e w _ 1 -s 

In (5) In (s) 


(9.3) 


The tables of Laplace transforms can also assist in evaluating the integral on the right hand 
side of equation (8. 3. 2.4). That is, if q(x)c an be written as ^(jc) = —a: h- r(jc) then we may write the 
transform as 

j 00 I 00 co 00 

/J J/c(jc) odf x) f{t) dx 1 = F(s)J k(x)s q ^dx = F(s)J ^ _ F(^)J ^~ vln ^^ r ^^ ln ^^(A)}dv.(9.4) 
[0 Jo 0 0 

Then, provided {e r ^ ln ^ 't«i is bounded and is of exponential order, its Laplace transform may 
be determined from the tables and the Laplace variable replaced by ln(.v). 
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Example application 

This example demonstrates the application of 
distributed order operators to obtain the analytical 
solution of the behavior of a viscoelastic material 
with spatially varying properties. We consider a 
plate constructed of thermorheologically complex 
material subjected to a temperature gradient T( x) 
and to a force F(t). Thus, the order q varies with 
temperature and therefore with the spatial 
dimension x. A free body diagram showing the 
distributed parameters is presented in figure 9-1. 
For any constant temperature element, the 
elemental stress- strain relationship [Hilfer p.35 1 ] 
is given by 

o(t)=Er\d?e(t), (9.5) 


/■« 

M') 

x 

1 1 



/=('> 1 

: 1 

+ 







- 1 r- 

+F(t) 


Figure 9- 1 Free body diagram 
for viscoelastic plate 


where a(t) is the applied stress and £*(/) is the resulting strain, q is the fractional order, with 
0<q< 1 , and E and r are material properties. Setting k q = Ez q we write for each incremental 
element and for a unit area 


*/(')— oD^Fit). 


(9.6) 


Then the total strain for the plate will be 


m /=i k <n 


(9.7) 


We will assume that the spatial temperature variation affects the order q such that q=c + dx , 
and in the limit as n -> «> , we have 




o‘.W 


oD t 


-(c+d .v) 


F(t)d x . 


(9.8) 


This is an initialized cumulative order distribution in terms of the independent variable x . For 
simplicity and clarity in this example, we will assume an initialization y = 0 , and / = 1 . Then 
taking the Laplace transform we have 


e(s) = j —U L{ 0 d; (c+d x) F(t)]d -v • ( 9 - 9 ) 

Again for simplicity we will assume for this example that E and r ? are unchanged by the 
temperature gradient, by taking k ? (,v)= 1 . Evaluating the transform gives 


i i 

f(^)= J , 5 _ ^ +rfv * F(s)d x= F(s). 

o o 


(9.10) 
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Performing the integration with F(t) taken as a unit step function, i.e., F(s) = l/s yields after 
some algebra 


e 



1 

s c+l ln(s) 


1 

c+l+d In (s) ' 


(9.11) 


The inverse transform may be found from the tables of transforms [Erdelyi, vol. 1, p.251] and the 
time domain solution is the sum of two Volterra type functions, thus 



t P+c 

r(p + e + l) 


dp - J 

c+rf+l 


j p+c+d 

r(/? + r + </ + l) 


dp\. 


(9.12) 


10. Operators of Multivariable Order and Distributed Variable Order 


10.1 Multivariable Order Operators 

In the past sections we have seen the order of fractional order operators vary with time and 
distributed over a (spatial) dimension x . We will now assign meaning to operators of the form 
0 D r ^’v)/(f), with order a function of multiple variables and combine the variable order and order 
distribution operators. We start by limiting consideration to uninitialized operators for clarity of 
discussion. Thus we define 


o d, | - r ^ v )) f( T ) dT > <?>0, (10.1.1) 

where the integral is to be evaluated with the variable .v treated as a constant. This definition 
would be used in applications where the order is an independent function of the (spatial) x 
dimension. An example might be to describe the local properties of a plate composed of different 
viscoelastic materials stacked along the * dimension, where each material has different (spatial) 
order properties. This definition is readily extended to multiple independent dimensions as 


/-<?( 


t, ' V; -' VJ /(0= j 


ri-if 


.V, ,.V; 


x n 


r(<?(- 


Xi , X~> JV 1 


„)) 


f(-r)dT. 


q> 0, (10.1.2) 


with each dimension treated as a constant. Application here generalizes the example above. That 
is to describe the local material properties for a viscoelastic solid composed of a three 
dimensional lay-up of thermorheologically complex materials and subjected to a temperature 
gradient T(.v). Thus, the order q varies with temperature and therefore with the x , y,c spatial 


dimensions along which each material has different order properties. 


For the above definitions the fractionally integrated function f(t) may also be generalized in 

an obvious way to be a function of x x ,x 2 x n ,t by use of the concept of the partial fractional 

integral. In this concept the function variables x l% x 2 ,...x n are treated as a constants in 
/(*!* *2 *••**#! >*)• Th e interested reader is referenced to Samko [1987 p.459] for discussion of the 
fixed order partial fractional integral. 


The (time) variable order fractional integral has been defined as (equation (4.1.2)) 




r( 9 M) 


-f(r)dT, 


q> 0. (10.1.3) 
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where for the “tailored” case we take q(t.z)= q(a t +bz) with a and b set by the application 
requirements. 

We now combine the t , (time) and x , (space) variable order operators and define 

( t __W( v.(f.r))-l 

/( ')-f lUCT /(t),,r - ,>0 ' (,01 ' 4) 

where again .vis taken as constant under the integration, and the notation (r,r) refers to the t , 
(time), variable order relationship as developed in section 4.0, e.g., (r.r)=> (at+bz). Application 
here might be to a viscoelastic solid composed of different complex thermo-viscoelastic materials 
along the .v dimension, where each material has different order properties, and those properties 
vary with a time varying temperature. For example let us assume the transient behavior relates as 
functions of ( t-z ), temperature is linear with time and baseline properties are linear (and 
independent) along the .v dimension then q(x,(t-z )) might be written as q(x,(t-z))= q(cx(t-z)). 
Extension of equation (10.1.4) to multiple dimensions yields 


r (t T Y^ V, * V: r)M 

0 v ' /('Wttt oT))) /(rVr - 

J rw(v,,.v : v„,(/,r))) 


(10.1.5) 

as the general uninitialized definition. The initialized definition is a direct application 

of property 

G2, thus for locally independent histories 





t > 0, q > 0, 

(10.1.6) 

where 



( .L L _ V/( V| ..v : v, f .(r ))- 1 

i/(/,v/(.Vi.v,,....v„,/)a,0,r)= J f(z)dz . 

J r(<y(.Vj , AS , (/, T ))) 

a 

r > 0,<? > 0. 

(10.1.7) 


The generalized variable order fractional derivative is defined as a generalization of equation 
6.2, thus we define 


o D, 


</(-> i 




/ (t H o D r" ( ' 1 ' " Vj 


t>0,q>0, (10.1.8) 


where q{ 1y,,.y 2 t)= m-u(x l ,.x 2 r), and m is taken as the least integer greater than the 

maximum value of q(x x ,x 2 /). 


10.2 Distributed Variable Order Operators 

The variable order and the order distribution operators may be combined to allow order 
distributions based on variable order operators. In a physical sense we are typically discussing 
time varying order distributions over some spatial extent. We limit discussion here to 
uninitialized operators for clarity. Then combining the results of equations (6.1) and (8.2.4) it is 
reasonable to write for the variable order direct approach cumulative distribution 

>’(') <?)= | *(<?) C d^ (,) y(l)dq , (l0.2.l) 

<?1 

or 

„., 3 tt£\k(q\ y(r\q)= } k(q) >'^ dZ d « ' ( 10 ' 2 ' 2 ) 

<?i c 
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where q(t)~ and the argument of u of the inner integral is selected as discussed in earlier 

sections. In similar manner the variable order independent variable (spatial) cumulative 
distribution is obtained by combining equations (6.1) and (8.3. 1.2) and is given as 

B 

f H v) c df xj] f{t)dx. (10.2.3) 

Alternatively using the form of equation (8.3. 1.3) we have 

nf'\k{q\ f(t\x)= j k(x) 0 df xj) f{j)dx, (10.2.4) 

0 

In terms of the Riemann-Liouville formulation this is written 

af !) (k(q\ f(t\ x)= j k(x) () dr J r))) f(r)dv dx. (10.2.5) 

where q(t,x) = m-u(xj) and the argument of u of the inner integral is selected as discussed in 

earlier sections. Then, for example using the Case 3 definition, for variable order fractional 
integration, i.e., q(t,r)= q(t - r), equation (10.13) becomes 


f k ( x ) / wrf(^)dr dx . (10.2.6) 

{ i r(«(.v, (/ - r ))) 


,(.t.(r-r)Vl 


In this case the Laplace transform may be written, using the result of equation (4.4. 4. 3), as 

Lh?< vj) (*(4 /(/U)^J/.-(.v).v m zj r ^ f *)) jf(5) dx. (10.2.7) 

where all initializations are assumed to be zero. Such formulations, as equations (10.2.5) or 
(10.2.7) may be applied to advantage. For example the viscoelastic problem associated with 
figure 9-1 could be subject to a time varying temperature, the form of equation (10.2.5) or its 
transform may be applied directly to determine the time varying cumulative strain for this 
complex system. To be determined for such a problem however, is which time variable order 
argument q{t,r) best describes the transient physics. 


The spatial strength spectrum may be time varying, i.e., k(x)~> k(x,t), also the integrand (input 
force) may be spatially distributed, i.e.,/(f)-4 f(t,x) then equation (10.2.5) generalizes to 

Qf x -' ) (k(xj). f(x,t\x)= j*(.r, /)„ d^'- x] f(x,t)dx . (10.2.8) 

0 

This form requires treating the fractional derivative (or fractional integral) as partial fractional 
derivative (or fractional integral); that is, taking x as a constant under the fractional operation. 

11. Summary 

This paper introduces the concepts of variable order fractional integration and differentiation 
where the order is a function of both time, / , and space, x . The variable order fractional integral 
is shown to allow considerable freedom in the argument of q , the order, which may be utilized to 
advantage in a variety of applications. The argument of q is also shown to control the manner in 
which the operator “remembers” the order history. Two types of memory measures have been 
shown. The first relates to the well understood, fading memory of the fractional operator, while 
the second type gives insight to the manner in which the variations in order within the operator 
are remembered. Measures have been introduced to quantify the memory effects. 
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Physical models for the three primary cases of variable integration have been shown. More 
significantly operational methods have been determined for the three primary cases and some 
applications of the methods have been shown. Extension from the primary cases to “tailored” 
operators was developed, and some properties of these tailored operators have been shown. 
Further research in this area is needed to derive a general operational method that applies to all 
the tailored variable order differintegrals. 

Initialization of the variable order operators has been shown. 

The fact that the variable order fractional operator is linear with respect to the input, /(/), and 
nonlinear with respect to q when order is considered as an input presents both problems and 
opportunity. The problem is the loss of the conventional composition law, which increases 
analytical difficulty. Until generalized composition laws are evolved for the variable order 
differintegral, the analyst is left primarily with numerical approaches. On the positive side, the 
dual linear / nonlinear nature of the operator may enable progress in the area of nonlinear 
systems. 

Two approaches to distributed order fractional operators have been presented. This distributed 
operator effectively combines all linear constant coefficient derivatives into a single operator. 
Laplace transform operations with distributed order operators have been derived and a 
viscoelastic strain application with a temperature gradient is shown. 

Finally, operators of multivariable order have been defined generalizing the variable order 
operators. Further, the effect of (time) variable fractional differentiation has been introduced into 
the order distribution, increasing its applicability. It is the author's belief that very many 
applications will be found for the new variable and distributed order operator definitions, 
forwarded here, in engineering and the sciences. 
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